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Abstract. This article deals with the numerical resolution of backward stochastic differ- 
| ential equations. Firstly, we consider a rather general case where the filtration is generated 

by a Brownian motion and a Poisson random measure. We provide a simulation algorithm 
based on iterative regressions on function bases, which coefficients are evaluated using 
Monte Carlo simulations. We state fully explicit error bounds. Secondly, restricting to 
the case of a Brownian filtration, we consider reflected BSDEs and adapt the previous al- 
gorithm to that situation. The complexity of the algorithm is very competitive and allows 
us to treat numerical results in dimension 10. 

Introduction 

> . 

Since the early nineties, there has been an increasing interest for BSDEs. Due to numerous 
■^j- \ contributions, we now have a better understanding of the features of these equations. This 

is still a very active field of research as this special volume may show. BSDEs have a wide 
\Q ■ range of applications, especially in finance, see for instance [8]. Our concern is rather 

related to the simulation of BSDEs, and regarding these aspects, the contributions in the 
literature are more recent and less definitive. Actually, the design of efficient algorithms 
that are able to solve general BSDEs in any reasonable dimension is far to be solved. 
For a review of different contributions, we refer to [9]. Here we present a new algorithm. It 
is designed for rather general BSDEs which are coupled with a forward Markov process X. 
The scheme that we propose has two ingredients. The first one is somewhat standard and 
consists in approximating the BSDE by a discrete-time backward dynamic programming 
equation. Then one needs to compute a sequence of regression functions: our choice is to 
project these functions on a finite-dimensional space spanned by basis functions, which 
coefficients are computed using a set of simulations of X. The details of the algorithm are 
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given in Section 1, in a generalized case where the filtration is generated by a Brownian 
motion and a Poisson random measure. In that case, X is a jump-diffusion process. We 
state explicit error estimates, with respect to the parameters of the method which are 
the time step h (used to derive the discrete-time dynamic programming equation; this 
parameter may be also used as a time step to simulate X), the number of basis functions, 
the number M of simulations of X, the threshold R needed to ensure a stability property 
(its impact can be usually neglected). 

We would like to emphasize several valuable features of our approach. Firstly, the explicit 
error bound allows us to choose optimally how the above parameters should vary together 
to achieve a given accuracy. Secondly, unlike the quantization approach [1], here the 
driver may depend on Z. Thirdly, we use only one set of M simulations to evaluate all 
the regression operators at once, while in [5] M x ^{discretization times} = 0{M/h) 
simulations are required. Finally, our scheme makes a very little use of the model for X. 
Indeed, on the one hand, for the implementation one just needs to know how to simulate 
approximative paths of X. On the other hand, for the analysis of regression errors, we 
use distribution- free techniques (see Gyorfi et al. [11]): hence a significant part of error 
estimates is model-free. We guess that this is a significant advantage, which allows us to 
consider general models (for instance, with the Malliavin calculus approach in [5], X is 
assumed to be essentially an elliptic diffusion). 

In Section 2, we consider reflected BSDEs, taking for X a diffusion process. We adapt 
the previous algorithm into three different variants, depending on how we handle the 
reflection. It alternatively gives lower and upper solutions. 

The scheme presented in this article has been first introduced and analyzed in [15]: therein, 
one deals with the framework given in Section 1. We refer the reader to the mentioned 
reference for the detailed proofs. Materials of Section 2 come from the PhD thesis of the 
second author [14]. 

1 Generalized backward stochastic differential equation 

We follow the presentation of Barles et al. [4]. Let (O, J 7 , (^)t,P) be a stochastic basis, 
where the filtration satisfies the usual conditions of right-continuity and completeness. 
We suppose that the filtration is generated by the two mutually independent processes: 
a IR 9 - valued Brownian motion W and a Poisson random measure /(/ on R + x E, where 
E = R'\{0} is equipped with its Borel field £, with compensator u(dt, de) = diA(de), such 
that {//([O, t] x A) = (n — ^)([0, t] x A)} t > is a martingale for all A £ £ with X(A) < +oo. 
A is assumed to be a cr-finite measure on (E,£) satisfying f E (l A |e| 2 )A(de) < +oo. We 
consider the Revalued jump-diffusion 

X t = x+ [ b{s,X s )ds+ [ a{s,X s )dW s + f f 0(s,X s -,e)fi(ds,de), (1) 

JO JO JO JE 

which is uniquely defined under the following assumption. 

(HI) The functions b(t,x) and a(t, x) are uniformly Lipschitz continuous with respect to 
{t,x) £ [0,T]xR d . For some constant c, the function (3 satisfies \(3(t, x, e)\ < c(l A |e|) 
and \(3(t,x,e)-(3(t',x',e)\ < c(|x-x'| + |t-t / |)(lA|e|) for any (t, x), (t, x') G [0,T]xI d 
and e G E. 
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Now consider the following generalized BSDE 



-dY t =f(t, X t , Y t , Z t )dt - Z t dWt - dL t , 



Y T = 4>(X T ), 



(2) 



where L is cadlag martingale orthogonal to W (with Lq = 0). Provided that the driver 
/ is Lipschitz continuous, there is a unique solution (Y, Z, L) (in appropriate L2-spaces, 
see [4] and [7] for details). To derive explicit error estimates, we impose to the terminal 
function <\> to be Lipschitz continuous as well. These two regularity assumptions are stated 
as follows: 

(H2) For any (t 1 ,x 1 ,y l , Zl ), {t 2 ,x 2 ,y 2 , z 2 ) G [0, T] x R d x R x R«, one has \f(t 2 ,x 2 ,y 2 ,z 2 )- 
f(h,xi,yi,zi)\ < C f (\t 2 - til 1 ^ + \ x 2 ~ xi\ + \y 2 - yi\ + \z 2 - zi\). The terminal 
condition <f> is Lipschitz continuous. 

1.1 Dynamic programming equation 

We first consider a time discretization of the equation (2). We denote the time step by 
h = and (tk = kh)o<k<N stand for the discretization times. For an arbitrary process 
U, set AUk = Ut k+1 — Ut k . One needs to approximate X by a Markov chain X N which 
can be simulated (take e.g. the Euler scheme [13]). Whatever the scheme we use, we only 
require that it converges to X in L 2 : 

(H3) sup < fc <jvE|X tfc - X%\ 2 -> as N goes to infinity. 

The discrete-time counterpart of (2) is given by Yj£ = <j)(X^ N ) and 



Y t N k =VM N k J + hV t J(t k ,X^Y» +v Z»), hZ» =E tfc (y t f +i A^), (3) 



where E tfe stands for the conditional expectation with respect to Tt k and * for the trans- 
pose. This is obtained by minimizing the difference E(Y t ^ +i + hE tk f(tk, X^, Yj£ , Z) — 
Y — ZAWk) 2 over J-% k -measurable squared integrable random variables (Y,Z). Our first 
result is the convergence of (Y N , Z N ) towards (Y, Z) in the BSDE -norm, as iV goes to oo. 
From Y N and Z N , one could also easily define a process L N and prove its convergence to 
L. The result below seems to be new in this general setting. 

Theorem 1 Under (H1-H2-H3), define the error 



where Y and Z are given by (3). Then, e(N) converges to as N —> oo. Furthermore, 
in the case of Brownian filtration ((3 = and L = 0) and when X N is the Euler scheme 
ofX, one has e(N) = O^ 1 ). 

Proof. In the sequel, C denotes a constant which value changes throughout computa- 
tions, but remains independent on TV. 

We begin by proving the result for maxo<fc<Ar E|y ifc — Y t ^\ 2 . Firstly, we know [7] that the 
solution (Y, Z, L) satisfies 




JV-l 



■t, 
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Then, from (2-3) we get Y tk - Y£ = E tk (Y tk+1 - Y* +1 ) + E tk X„ Y„ Z s ) - 

f(t k , Xg, Y£ i , Z^)}ds. A combination of Young's inequality (a + b) 2 < (1 + 7 /i)a 2 + (1 + 
^-Jb 2 (with a parameter 7 > to be chosen later) and of the Lipschitz property of / gives 

1 ftk + l 

E\Y tk - Yt N k \ 2 <(1 + jh)E\E tk (Y tk+1 - Y t N k+i )\ 2 + C(h + -)E / \Z S - Z t N fds 

+ C(h+-)(h 2 + E\X s -X t N k \ 2 ds+ E\Y s -Y£ +1 \ 2 ds). (5) 

If J tk Jtk 

Now define Z tfc by /iZ tfc := E tfc Z s d S = E tfc ({Y tfc+1 + f(s, X s , Y s , Z s )ds}AW*). 
Clearly 



E \Z s -Z t N J 2 ds = E \Z s -Z tk \ 2 ds + hE\Z tk -Z t N k \ 2 . (6) 

Jt k Jt k 



The Cauchy-Schwarz inequality yields \E tk {{Y tk+1 - Y^ +i }AW l>k )\ 2 < h{E tk (\Y tk+1 - 
Yt N k+1 ?) ~ \V tk (Yt k+1 ~ Yt N k+1 )?} and consequently 

hE\Z tk - Z t N J 2 <CE{E tk (\Y tk+1 - Y t N k J 2 ) - E\E tk {Y tk+1 - Y t N k+i )\ 2 } 

ftk + l 

+ ChE f(s,X s ,Y s ,Z s ) 2 ds. (7) 

Jt k 

Plugging (6-7) into (5), we get: 

1 ftk+l 

E\Y tk - Y t N k \ 2 <(1 + jh)E\E tk (Y tk+1 - Y t N k J\ 2 + C(h + -)E / \Z S - Z tk \ 2 ds 

7 Jt k 

+ c(h + -)(h 2 + / E\x s -x t N j 2 ds+ / E|y s -y,f +i | 2 d s ) 
7 Jt k Jtk 

+ C(h + l -)E{E tk (\Y tk+1 - Y t N k J) - \E tk (Y tk+1 - Y t N k J\ 2 } 
1 ftk+i 

+ Ch(h + -)E / f(s, X s , Y s , Z s ) 2 ds. 
7 Jtk 

Now write E| Y-Y t f + J 2 < 2E\Y S -Y tk+1 \ 2 +2E\Y tk+l -Y t N k J 2 and analogously for X s -X? k , 
take 7 = C: for h small enough, it gives 

E\Y tk - Y t N k | 2 <(1 + Ch)E\Y tk+1 - Y t N k J 2 + Ch 2 + Ch m^E\X tk - X? k \ 2 



+ CE / \Z S - Z tk \ 2 ds + C E\X S - X tk \ 2 ds 

Jtk Jtk 
rtk+i ftk+i 
+ C E\Y S - Y tk+1 \ 2 ds + ChE / f(s, X a , Y s , Z s ) 2 ds 

Jtk Jtk 

r t k 



and by Gronwall's lemma maxo<fc<7v E\ Yt k — Y t ^\ 2 < Ch + C maxo<fe<7v E\ Xt k — 
Xt N J 2 + CJ2^Ef t [^{\Z s - Z~J + \X S - X tk \ 2 + \Y S - Y tk+1 \ 2 }ds. The contri- 
bution J2k=o It k +1 ^\Ys — Y tk+1 \ 2 ds is a 0(h): indeed it is upper bounded by 
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3 E^To 1 ItT 1 Mtk+i - s) f l s k+1 Ef(u, X u , Y u , Z u ) 2 du + 3 Y.k=o J? +1 ^ nZ u \ 2 du + 
3J2k=o ftk +1 dsE([L] tk+1 ~~ Ms)) which equals a 0(/i) owing to the a priori estimates (4) 
on (Z,L). In the same way, the contribution related to X s — X tk is of order 0(h). Finally, 
it gives 

N-l rtk+i 

max E\Y t . - Y, N \ 2 < Ch + C max E|X t . - X f N \ 2 + C V E / |Z S - Zj 2 ds. 

0<k<N fc fc 0<k<N fc ^ L 

— k=0 Jtk 

Without extra assumptions, the above approximation's error related to the predictable 
process Z converges to 0: combining this with (H3), we conclude max <fc<Ar E|Yi fe — 

In the Brownian filtration case ((3 = 0) and when X N is the Euler scheme of X, clearly 
E|X tfe — Xj^\ 2 = 0(h) uniformly in k. Furthermore, Zhang [18] establishes that the error 
on Z equals 0(h). Hence, max <fc<Ar E\Y tk — Y t ^\ 2 = 0(h). 

Using the same techniques, it is easy to derive an estimate for E Efc^o 1 ff k+1 \^tl ~ ^t\ 2< it 
from that on Y — Y (see also the proof of Theorem 2 in [9]). We omit further details. □ 



1.2 Description of the algorithm 

At some point in our approach, the numerical solution of (Y, Z) needs to be upper bounded, 
especially to analyze the empirical regression errors. This is a theoretical reason but we 
also mention that it seems to have some importance in our numerical experiments. To 
get a bounded solution, we consider a threshold R = (R , R u . . . , R d ) € (a priori 

with large coordinates) and this threshold defines 

[A^,i = ( - R Vh ) V AW,,* A (R Vh), 
f R (t,x,y,z) = f(t,-R! Vxi Ai?i,--- ,-R d V x d AR d ,y,z), 
<p R (x) = <P(-R l V xi A Ri, ■ ■ ■ , -R d \/x d A R d ), 

M v (x) = -C y (R) V ^(x) A C y (R), &} z (x) = V ^ (x) A 

The constant C y (R) is explicit and defined later (see equation (9) below). 
It is easy to see that the conditional expectations in (3) are regression functions of Xj^; 
hence, it is natural to approximate them by their projections on finite-dimensional bases. 
For the scalar component of Y and the q components of Z, at each discretization time 
we take q + 1 deterministic function bases {Pl,k(')) o<i< q - ^ or sa ^ e °f convenience, the basis 
Pl^k of size Ki k is considered as a iQ^-dimensional" vector of scalar functions. A function 
F in the vector space spanned by the basis p\^ is given by its coefficients a and one uses 
the short notation F(-) = a ■ pi,k(-)- 

We denote by (X tk ' m )i< m < M ,o<k<N and (AW™)i< m <M,o<fc<iV-i the M independent sim- 
ulations of (X^) < k < N and (AW fc )o<fc<iV-i- 

With these notations, we are in a position to define yj?' R ' M (X^) and zfy R ' M (X%) (the 

approximations of Y£ and Zg fc ), where the functions ' R > M (■) and zfy R ' M (•) are given 
as follows. 
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Initialization : for k = N set y^ ,R,M (■) = 4> R {-)- 

Iteration : for k = N — 1, ■ ■ ■ ,0, solve the q least-squares problems 1 : 



M • r 1 \ ~* i N,R,M i v N,m\ 

a l:k = arg mf — ^ | y k ^ (X tk ' +i ) a ■ Pl , k {X tk ) | . 



M ^ ,ak+L v tk+1 ' h 

m=l 



Put 



^ iN fc R ' M (') = i a ik ' Pi,k(-)]z • Then compute aff k as the minimizer of 



1 M 

1 l N,R,M, v N,m\ , ir-Ru v N,m N,R,M / v N,m\ N,R,M tv N,m\\ (vN,m\\2 

jjjZJtffc+i ( X t k+1 ) + h f (tk,X tk ' ,y k ^ {X tk ' +i ),z hk - {X th > ))-a- Po , k (X tk > - 



m=l 



over a. In the above definition, we write f R (t k ,x,y,zi) for f R (t k , x, y, (z;)i</< 



Put 



N,R,M i \ r M I m 



1.3 Convergence analysis 

1.3.1 Error from the threshold R 

We first compare the approximative solution (Y N ,Z N ) defined by (3) with (Y N,R , Z N,R ) 
defined by: 

Ytf* = Efc(y£f) + hK t J R (t k ,Xi[,Y t ^ R ,Z^ R ), hZ^ R = E tk (Y» ; R [AW k ]*J, (8) 

and Y t * >R = <p R (X t N N ). In [15], we prove that (Y N ' R , Z N > R ) converges in L 2 to (Y N , Z N ) as 
R goes to infinity. A rate of convergence is also available under the additional assumption 
sup <fc<JvE|X^|P < C p (l + \x\p) for some p > 2. Namely, provided that R { = iV 2 /^ 2 ) 
{i = l,--- , d) and i?o = cy / log(AT) (for c large enough), one has max <fc<jv E\Y t ^ ,R — 
Y*\ 2 + hEJ2k=o \ z t k ' R ~ Z t k ? = OiN- 1 ) as in Theorem 2. Hence, especially if p is not 
small, reasonable values of R yield a very good approximation. 

The new dynamic programming equation (8) has the main advantage to give bounded 
solutions. Namely one has \Y t f R \ < C y (R) and Vh\z£' R \ < C y {R) where 

C y (R) = e^ +1± f^ {sup |<^(x)| 2 + 2Ti±^ sup \f R (t, x, 0, 0)| 2 ) (9) 

I x 7 t,x J 

and 7* = AqC 2 (see [15] and [14]). Note that the Loo-norms related to <fi N ' R and f R can 
be easily computed on usual examples. 

1.3.2 Error from the projections on bases and the simulations 

The following extra assumption (H4) is sufficient to prove that the functions y k ' R {) 
and y/hz?' R (-) (defined by y^' R (X t N k ) = Y t f R and z£' R (X t N k ) = Z% R ) are Lipschitz 
continuous, uniformly in h and R (see [15]). 



1 For the numerical resolution of least-squares problems, see [10]. 
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(H4) Denote by X^ k ' k °' x the Markov chain Xj^ starting at x at time t ko . There is a 
constant C > such that 

a) E\X^ k °> x - X?/°' x '\ 2 + nX^ o k + °; x - X%' k °; x '\ 2 < C\x - x'\ 2 for any x and x\ 
uniformly in ko and N. 

b) E\X?' k °> x - x\ 2 < Ch{l + |x| 2 ) for any x, uniformly in k and N. 

This assumption on X N is not restrictive since this property holds for X under (HI) (and 
usually, for the Euler scheme as well). The next result states a general upper bound for the 
error between (Y t ^ ,R , Z^' R ) and the solutions (y^' R,M (Xj^), z^ ,R ' M (Xj^)) computed with 
our scheme, in terms of the basis functions and the number of simulations. Specifications 
of the upper bounds are given later, according to the choice of bases. 

Theorem 2 (see [15]) Assume (H1-H2-H3-H4) and denote by Kf\. the random variable 
given by the rank of the matrix of size x M which columns are (pi^(xj^' m )) m (note 
that KM < Kik). By convention, we set Kq^ = 0. 
Then, there exists a constant C such that for any (5 G]0, 1] one has 



N-l 

%N,R, V N\ N,R,M / vJV\i2 , uv? I N,R, V N\ N,R,M /V 
, ,. v k ( x t k )-y k i x t k )\ + ^ E 2JV ( x t k )~ z k ( x . 

— fc=0 



M 

k=0 1=0 



N-l 



+ c zZ {^nyk' R (x t N k ) - *-PoA x t N k )\ 2 + Y,^n^ R (x t N k ) - a . Pl , k (x t N k )\ 2 } 

k=0 1=1 

, r Cy(Rl N ^( w , M Mh^ 2 C CyimK^Yl 

+ c—^ g |e(^ *M- nCy{R)2K M ) *MCK , k+1 log w )) 



+ hRlKfi exp( 9 9 M ) exp(CK k+1 log yK ' + , l > kJ ) 

V i,k 72C y {R) 2 R 2 ) K l M k J v ' + " 

CCy(R)^ , Mh^ 2 



+ exp(CK ,fc log j+2—) exp( 



h ep ,~*k nCy {R) 2 ' 

The parameter (3 should be chosen optimally to get the sharpest upper bound, see Para- 
graph 1.4. Theorem 2 improves our previous results [9], where the statistical errors are 
estimated in terms of the fourth moments of the L2-orthonormalized basis functions. These 
quantities are generally unknown, which makes the adjustment of the bases, h and M dif- 
ficult. In order to interpret terms of the upper bound above, we first recall a well-known 
result in non-parametric regressions. 

Theorem 3 (see [11] Theorem 11.1 p. 184)- Let (U,V) be two random variables taking 
values in $l d x 3? and let F be a K -dimensional linear space of functions. Consider the 
problem of estimating the regression function v(u) = E(V\U = u), using M independent 
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realizations (U m , V^,)i< m <M of (U,V). Approximate v(-) by its best approximation in F 
using an empirical L,2-norm related to the M observations: 

1 M 

v M {-) = arg mf — ^ \V m - f{U m )\ 2 . 
m=l 

Under the assumption a 2 = sup u Var(V|f7 = u) < oo, one has 

f 1 M \ K 

E ( M E \< U m) ~ VM(U m )\ 2 \ < a 2 - + minE|/([/) - v(U)\ 2 . (10) 

^ m=l ' 

This is the prototype of model-free results since one only requires the conditional variance 
to be uniformly bounded. 

In our case, at each tk one solves a regression problem, which yields errors of type (10) 
(compare with the first and second lines of the r.h.s. of Theorem 2). The boundedness 
of (Y N ' R , \fh~Z N ' R ) ensures that the relative conditional variances are uniformly bounded. 
The additional factor log(M) in the first line occurs because one estimates R\y%' R (X%) - 

N,R,M / v"JV\ |2 • +. i r 1 I N.R, v N,m\ N,R,M , v N,m\\2 ■ n n\ 

Vk ( X tk)\ instead of E s £ m=1 \y k [X tk ' )-y k ' (X t ^ )\ 2 as in (10). 
However, the scheme to currently analyze is much more delicate than in Theorem 3, 
because the regression problems are all correlated: indeed, V m should be given by 
Uk+i' M '(Xtl'™) which is not independent of V m i (m / m!) because of the random function 

yk+i' M (')■ ^ ne trick consists in covering the class of functions where y^^ M {-) is, by a 
finite number of balls centered at deterministic functions. The radius of the above balls 
depends on h and on an extra parameter (5. The covering number is estimated using the 
Vapnik-Chervonenkis dimension. All these extra contributions are the other terms in the 
r.h.s. of Theorem 2. 

The analysis of the algorithm's error turns out to be quite intricate and actually, we guess 
that our estimates of Theorem 2 are not optimal. However, we do not have good ideas 
to improve them. Before performing an analysis of complexity in Paragraph 1.4, we give 
a modified algorithm which complexity is the same and which theoretical upper bound is 
better (paradoxically, both schemes behave similarly on numerical tests). 



1.3.3 A modified algorithm 

As mentioned before, some terms in the upper bound of Theorem 2 come from the lack of 
independence between least-squares problems at different discretization times tk- These 
terms are those from the third and fourth lines in the r.h.s. of the upper bound. We 
present below a variant of the algorithm of Paragraph 1.2, which removes these two terms 
(for details, see [14]). 

— ► Initialization : for k = N set y^' R,M (■) = 4> R {-)- 

— ► Let k < N— 1. For each path m, simulate (X^'™, AW 7 "™) which are, conditionally to 
Xj^ ,m , an independent copy of (X^™, AWJJ 1 ) (and independent of everything else). 
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The coefficients af^ and a^ k are now computed as 

a, )fc = argmf— 2Ji/k+i' (^t k ' +1 ) £ «-Pi,fc(^t fc ' )h 

m=l 
M 

a o,k ~ ar g!ni M 2^ Wk+l y^t k+1 ) 

m=l 

, u*Ru v N,m N,R,M , -$-N,m\ N,R,M , v N,m\\ , v N,m\\2 

+ hf (t k ,X tk > ,y k ^' (X tk ' +i ),z lk > (X tk ' )) - a ■ p , k (X tk ' )| . 
As before, the functions y k ' R,M (■) and zfy R ' M (•) are given by [a^, -po.fcOlj/ an d 

[°S •«,*(■)].■ 

Note that this algorithm requires to draw only twice more simulations. 
1.4 Accuracy and complexity 

Now we focus on the error between the discrete-time dynamic programming equation 
(3) and the solutions computed by our initial algorithm, its modified version and also 
the Bouchard- Touzi's algorithm [5]. Our aim is at discussing the trade-off between the 
accuracy and the complexity (i.e. computational time). 

1.4.1 Initial algorithm (see Paragraph 1.2) 

From Theorem 2, we can derive how to make N, M and the number of basis functions 
vary together. In this discussion, we neglect the influence of the localization parameter 
R, which is supposed to be large enough from the beginning (see paragraph 1.3.1). As 
already observed in [9], local basis functions take advantage of the Lipschitz property 
of functions y k ' R (-) and Vhzfy R (-) (see Paragraph 1.3.2). Let us consider the simplest 
example of a local function basis, i.e. the hypercubes basis, already used in [9] and 
still denoted HC here. To simplify, pi^ does not depend on I or k and its size equals 
K. Choose a domain D C M. d centered on x, that is D = fjf=i ] x * ~~ a i x i + °]> an d 
partition it into small hypercubes of edge 5. Thus, D = Uj 1; ... ,i d Di lt ... j d where Di lt ... t i d = 
]x\ — a + i\5 : x\ — a + {i\ + 1)5] x • • • x]xd — a + id5, Xd~ a + (id + 1)<5]- Then we define p\^ as 
the indicator functions associated to this set of hypercubes: Pi : k{~) = (l-D n ... i d ('))i 1 j d - 
With this particular choice of function bases, we can make the projection error of Theorem 
2 explicit and refer to [9] for details: 

miK(\y^ R (X t N k )-a-p , k (X t N k )\ 2 ) < C{6 2 + C y (R) 2 F(X t N k eD c )}. 

As for the impact of the threshold R, F(A^ £ D c ) becomes negligible with respect to 
the other errors if we choose D big enough (a feature which is confirmed by the numer- 
ical experiments). Thus and as far the projection errors are concerned, to get a global 
(squared) error of order ft' 3 we have to choose 5 « h~z~ , or equivalently a number of 

d(/3+l) 

basis functions K « ft 2 (considering a fixed domain D). Regarding now the num- 
ber of simulations M, to avoid an explosive upper bound in Theorem 2, one should take 
M « (7/ l - rf (/ 3 + 1 )-(/ 3 + 2 ) log (ft, 4 5~) for a constant C large enough (here, the ranks 
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Kfj, have simply been upper bounded by K). 

The dominant term of the algorithm's complexity C associated to this choice of function 
basis is C = NMdlog(K), which corresponds to determine in which cells the simulation 
fall (this is the cost of a nearest neighbor algorithm in a tensored grid, i.e. 0{d\og(K)) 
for one path at a given time). Hence, up to logarithmic factors, the complexity equals 
C = 0(/ l -i-^(/3+i)-(/3+2)) ; while the square d error is of order hP = 0(C-/ 3 /( 2 +(/ 3 + 1 )( fl! + 1 ))). 

The optimal value of (3 G]0, 1] is (3 = 1, for which the squared error is of order 
h = 0{C- l /( 2d+ ^), while 5 « h and M « Ch~ 2d - 3 log(l//i). 

1.4.2 Modified algorithm (see Paragraph 1.3.3) 

The complexity is unchanged. But the error bounds are better since the terms from the 
third and fourth lines in the r.h.s. of Theorem 2 have disappeared. Hence the asymptotics 
on M is less stringent: indeed, we get M « Ch~ d ~ 3 log(l//i), 5 ~ h in order to achieve a 
squared error of order h = 0(C-V(*^)). 

1.4.3 Bouchard- Touzi's algorithm [5] 

We compute the complexity of their algorithm in the most favorable case where X is 
a geometric Brownian motion. Otherwise in a more general diffusion framework, the 
algorithm is heavy to implement and its complexity more difficult to evaluate because of 
the necessary calculation of Skorohod's integrals. 

In this algorithm, one needs N independent sets (Mk)i<k<N of M simulated paths of X N 
(one set for each discretization time). At each discretization time and for each path 
in the set Mk, a calculation involving the M paths of the set Mk+i is performed. This 
leads to a complexity C = 0(NM 2 ). The squared error associated to this complexity is 

2+ ^ 

given by Theorem 6.2 in [5] and is of order + N ^ ■ Expressing the squared error as a 

i 

function of the complexity, we find C 13 + d . 

1.4.4 Summary 



Initial Algorithm 
(Paragraph 1.2) 


Modified Algorithm 
(Paragraph 1.3.3) 


Bouchard-Touzi's Algorithm 
[5] 


i 

C 4+2d 




i 

(if A=Brownian motion 
or geometric BM) 



Table 1: Squared error with respect to the complexity C. 



It turns out that in the geometric Brownian case, our algorithm is more efficient than 
Bouchard-Touzi's Algorithm for d < 9 and less efficient otherwise. But the geometric 
Brownian framework is very favorable for a Malliavin calculus approach, since the integra- 
tion by parts formulas are especially simple and this is no more true for general diffusion 
models. With our approach, the complexity is independent of the model. 
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Furthermore, the modified algorithm is very competitive in any case, although it is less 
natural (because of the extra simulations). We have observed in practice that the initial 
and modified algorithm differ very little from each other. This suggests that the upper 
bound in Theorem 2 may be not optimal (and the resulting complexity as well). We hope 
that we could address these issues in future works. 

2 Reflected backward stochastic differential equation 

In this section, we restrict to a Brownian BSDE reflected on one barrier. We show that 
the results presented in the non-reflected case can apply to the reflected case. Actually 
we present three approximation methods, which we call in short max method, penalization 
method and regularization method. 

The framework is analogous to Section 1, without the Poisson random measure process. 
Namely, we suppose that the forward process is a diffusion SDE X t = x + £ b{s, X s )ds + 
Jq a(s, X s )dW s under the following assumption 

(HI') The functions b and a are continuously differentiable with uniformly bounded deriva- 
tives w.r.t. x. The matrix aa* satisfies the ellipticity condition go* > eld with e > 0. 

Regarding the reflected BSDE (RBSDE in short), we consider a continuous function <p 
which defines a continuous adapted obstacle 4>{t, X t ). The RSBDE (Y, Z, K) is solution of 

' Y t = 4>(T, X T ) + ff f(s, X s , Y s , Z s )ds + K T -K t - Z s dW s , 
< Y t >cf>(t,X t ), (11) 
K is continuous, increasing, Kq = and f Q 7 (Y t — <p(t,X t ))dK t = 0. 

To ensure existence/uniqueness of the process above [6] and further approximation results, 
we impose stronger conditions on the driver and the obstacle: 

(H2') For any (tj, x u y u z ± ), (t 2 , x 2 , y 2 , z 2 ) £ [0, T] x R d x R x R q , one has \f(t 2 , x 2 , y 2 , z 2 ) - 
/(ti,xi,yi,zi)| < C f (\t 2 - t^ 1 / 2 + \x 2 - xi\ + \y 2 - yi\ + \z 2 - zi|) and \((>(t 2 ,x 2 ) - 
0(*i,asi)| < Ct(\h ~ h\ 1/2 + \x 2 - xi|). 

When the driver / is linear, i.e. f(t,x,y,z) = —rty — z Q% where 9 is the market risk- 
premium in finance, it is well known that Yt equals the price of the American option with 
payoff P t = <f>(t,X t ): Y t = ess sup stopping timcs r6 [ t)T ] E Q(e _ & rsds P T \^t)- 

2.1 Approximation procedures 

We denote by X N the Euler scheme of X, based on the time net (tk = kh)o<k<N where 
h = jj. As before, one needs M independent simulations of (X N , AW). We propose three 
procedures to approximate (Y, Z), solution of the reflected BSDE above. 

2.1.1 Max method 

This method is standard and simply consists in taking at each time tk the maximum 
between the obstacle and the expected value of the non reflected BSDE. The resulting 
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approximation (Y N , Z N ) is defined by: 

[Y»= cP(t N ,X? N ), 
Y t N k = max (<f>(t k , X t N k ), E tk (Y* +1 + hf(t k ,X t N k , Y t N k+i , Z t N k ))) , < k < N - 1, 
K= ^t k (Y t ^ +1 AW*),0<k<N-l. 

Ma and Zhang [16] prove that the rate of convergence of such approximation is A -1 / 4 
when the obstacle function is smooth (i.e. of class C 1,2 ). In [14] we extend this analysis 
to less smooth obstacles in order to handle for instance the case of American put contract 
(i.e. <f)(t,x) = (K — e x ) + if X is the log-price process). The relevant assumption of 4> is 
the following one: 

(H3') The function (ft is of class C 1 w.r.t. the time variable and uniformly Lipschitz contin- 
uous w.r.t. the space variable. Moreover, <fi(t, X t ) and <fi(t, X^) satisfy the following 
ltd expansions: 

4>(t, X t ) =0(0, x) + f U s ds + f V s dW s + A t , 
J J 

<f>(t, X t N ) =0(0, x) + f U^ds + f V s N dW s + A?, 
Jo Jo 

where A, A N are increasing, continuous, integrable and such that the measures dAt 
and dA^ are singular w.r.t. dt and where U, U N , V and V N satisfy the follow- 
ing integrability condition sup t<T E(|C/ t | p ) + sup t<T N E(|f// V | p ) + sup t<T E(| Vt\ p ) + 
sup t < TiJV E(|V/Y) < oo for p >~2. 
Theorem 4 (see [14]) Under Assumption (Hl'-H2'-H3'), one has 



N-l ,.t k+1 

ma^E|y t f - Y tk \ 2 + £ j E\Z t N k - Z t \ 2 dt < (7(1 + M 4 )A~ 



1/2 



0<fc<iV it, 

k=0 Jtk 



Note that the convergence rate is lower than in the non-reflected case (N 4 instead of 
N~?). The algorithm to approximate (Y N ,Z N ) is very similar to the one of Paragraph 
1.2 (with analogous threshold functions). 

-. N,R,M i \ iR/4. \ 

2. At a given discretization time t k (0 < k < N — 1), the coefficients af^ (1 < I < q) 

[AWj ' 



1 M [&W m ] w 
are defined as the minimizers over a\ of — ^ \y k ^ ,M (X^'™) ai.p™ k 2 

m=l 

Then put z^ M (.) = [aff k . Pl , k U-)- 
3. a^f k is finally defined as the minimizer over q of 
1 M 

V* I N,R,M fv N,m\ . utRn vN,m N,R,M / v N,m\ N,R,M , v N,m\\ rn 

Jj2^\y k +i ( X t k+1 ) + hf {t k ,X tk > ,y k ^ {X tk ' +i ),z k ' {X tk > ))-a .p 0>k 



m=l 

Then put y k ' R ' M (■) = max(^(t fc , •), [< fc .Po,fc]v(-))- 

The error between (Y t ' ,Z t ' ) and (y k ' ' (X^),z k ' ' (A^)) can be analyzed exactly 
as in Theorem 2 and actually, we get the same estimates. 
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2.1.2 Penalization method 



This method is based on the ideas developed in [6] to show the existence and uniqueness 
of the solution to (11): the reflection is handled via a penalization term which is very large 
only when Y is below to the obstacle. Namely, the solution (Y, Z, K) is approximated by 
(Y n , Z n , K n ) where (Y n , Z n ) is the solution of the standard BSDE 

Y t n = cf>(T,X T ) + [ T f(s,X s ,Y:\Z:)ds + n C (Y? - </>(s, X a ))-ds - C Z?dW s , 
Jt Jt Jt 

and K n is defined by = n J Q * (Y s n — (f>(s, A s ))_ds. It is thus shown in [6] that 
(Y n ,Z n ,K n ) tends to (Y, Z, K) as n tends to infinity with the property Y n | Y. In 
general, one does not know the convergence rate (some results are available when / does 
not depend on Z, see [12]). At a fixed n, (Y n , Z n ) is the solution of a standard BSDE and 
can thus be approximated by the algorithm proposed in the first section: 

1. y% N ' R ' M (-) = <i> R (tN,-) 

2. At a given discretization time tk (0 < k < N — 1), the coefficients (aff,) are de- 

n , ... r- r 1 i n,N,R,M , v N,m\ 

hned as the mimmizers of ini — Wk+i \ t k +i > J. a l-Pl,k\ ■ Ihen put 

m=l 

<f AM (-) = K-piM-)- 

3. a^ k is finally defined as the minimizer of: 

1 M 

■ r 1 I n,N,R,M i v N,m\ , , rR u v N,m n,N,R,M , v N,m\ n,N,R,M / v N,m\\ 

^mZJ^+i ( X t k ' +1 ) + hf (t k ,X tk > ,y k ' +1 ' (X tk ' +i ), Zl ; k ' (X tk ' )) 

m=l 

■ 7 / n,N,R,M i v N,m\ 1R14. vN,m\\ m |2 

+ "Mj/fc+i (tk,X tk > ))_ -a .p ,fel • 

Then put yl' N ' R ' M (-) = [a$ k .po,k] y (-)- 
2.1.3 Regular izat ion method 

The last idea is based on the results of [2], which exploit the specific form of the in- 
creasing process The authors consider the case where the obstacle 4>(t,Xt) is a 
cadlag semimartingale satisfying (p(t,Xt) = (f>(0,x) + Jq U s ds + J^V s dW s + At where 
E J Q {\Vt\ 2 + |£/t| 2 }eZ£ < 00 and where A is a cadlag adapted increasing process with 
E|Ar| < 00. In this case, a solution (Y, Z, K) to (11) is also a solution to 

'Y t = <f,{t,X t ) +/f {f(s,X s ,Y s ,Z s ) +a s l Ys=4l{S;Xs) [f(s,X s ,4>(s,X s ),V s ) + U s ]-}ds 

-J t T Z s dW s , 
Y t > <f>(t,X t ) 

with E J T {|Yf| 2 + \Zt\ 2 + \at\ 2 }dt < 00. To show the existence and uniqueness of the 
solution to the equation above, the authors use a regularization method based on C°° 
functions tp n (with < ip n < 1) satisfying ip n (x) = 1 if \x\ < ^ and ip n (x) = if \x\ > 
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which are mollifiers of (p(x) = l x =0- It turns out that the solution (Y n , Z n ) to the following 
BSDE converges to (Y, Z) as n goes to infinity with the property that Y n j Y: 

Y t n =4>(t,X t ) + £f( S ,X s ,Yp,Z^+Myr-Hs,X s ))[f(s,X s ^(s,X s ),V s ) + U s ^ds 

z n s m s . 

Note that it gives an upper approximation for Y while the two previous methods give lower 
approximations. We can thus hope to define an approximation procedure which gives an 
upper estimation for the Y (or equivalently for the prices of American options) . To keep 
this property, the time-discretization (Y N , Z N ) is performed in the following way: 

'Y»= Ht N ,X t N N ), 

^(Y^ + hfi^x^^z^ + hMY^ - 4>(t k+ i,x t N k+1 )) 

<t>(t k+1 ,x" )-6(t k ,Xf) \ 
[f(t k , X», cf>(t k+ n X t N k+1 ), V") + tfc Y A -) , 

E tk (Y t lAW*), 
[hV»= E tk (cf>(t k+1 ,X t N k+i )AW*). 

Then, the scheme of Paragraph 1.2 becomes 

-. n,N,R,M i \ iRu \ 

2. At a given discretization time t k (0 < k < N — 1), the coefficients (af|) are defined 

■ ■ ■ r 1 l n,N,R,M , v N,m\ 

W „ .9 _ 

the mimmizers over a, of — ^ \y k ' +1 ' ' { x t k ' +1 ) ^ a l-Pl,k\ ■ Then P ut 



as 



m=l 



M 



1 M 

3. Define 0% = urging- £ |^(t fc+ i, 



[AW, 



m=l 



h 



n m |2 



4. Finally a^ k minimize over «o 
1 M 

V Ln^,iJ,M, y JV,ms , , fj R/. yN,m n,N,R,M / yN,m\ n,N,R,M / yN,m\\ 



m=l 



+ hMv k 'T' M (Kp ~ ^ R (tk + i,x^))[f R (t k , x?> m , <t> R (t k+1 ,x»p, \pff k .p% k ] v ) 
*^x%™)-<t> R {t k ,x^) 



+ 



h 



Then put y k ' N ' R ' M {-) = [a$ k .po,k] y (-)- 
In practice, ip n is a piecewise linear function with ip n (x) = 1 if \x\ < ^, (f n (x) = if 
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2.2 Numerical results 



We present here numerical results for American options written on one, three and ten- 
dimensional assets. 



2.2.1 Dimension d = 1 

First, we take the same examples than in [5] where the authors consider the case of a usual 
Black-Scholes model in dimension d and the payoff 4>(x) = (K — (JliLi x i) 1 )+- 
We first take d = 1, the maturity T of the option is one year, the interest rate r = 0.05, 
the volatility a = 0.15, the strike K = 100 and Sq = 100. The reference price obtained 
by a PDE method is 4.23 according to [5] while the price obtained in [5] by a Malliavin 
calculus approach to solve the RBSDE is 4.21. For our algorithm, we use the basis HC 
of hypercubes with size 5 (see Paragraph 1.4.1). We test the max method, for different 
values of N and 5 and we let M increase. On Figure 1, we observe that the bias decreases 
when N and 5 increase, as predicted in Theorem 2. 




N=10,deltii=10 
N=20,delta=5 
N=20,delta=! 
N=50,delta=l 
Real Price 



1. 10 s 2. 10 s 3. 10 s 



N=20,de!ta=l 
N=50,de!ta=l 
Real Price 



Figure 1: American put with d = 1 



2.2.2 Dimension d = 3 

We deal with the same example with d = 3. We test the regularization method with 
N = 32, 5 = 9 and n = 2, the max method with N = 44, 5 = 7 and the penalization 
method with iV = 60, 5 = 2 and n = 2. For all of them, we use the basis HC. On Figure 
2, we observe that the three algorithms give good approximations of the American prices. 
Nevertheless, we mention that when parameters (N, 5, M) simultaneously change with the 
regularization or penalization parameter n, the methods may behave very differently and 
we refer to [14] for a more complete numerical study. 

2.2.3 Dimension d = 10 

Finally, we give an example with d = 10 = 2p taken from [3]. The model is still a Black- 
Scholes one : = (r — [Andt + aidW} for 1 < / < 10. The pay-off (or obstacle) is 
max(xi ■ ■ ■ x p — Xp+\ ■ ■ ■ X2 P , 0). The interest rate is r = 0, the dividend rates fi\ = —0.05 
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Regularization 

Penalization 

Max 

Reference price 



Figure 2: American Put with d = 3 
and [i\ = for I > 2, the volatility a\ = 77=, the maturity is T = 0.5 and the initial values 

2 ■ 2 

are Sq = 40^, 1 < 2 < p and <Sq = 36 3 for p + 1 < i < 2p. The reference price used in [3] 
is the one derived by [17], using a PDE method, and is equal to 4.896. 
We use the max method with a slightly different function basis : instead of just using 
indicator functions to approximate Y, we define also on each hypercube polynomials of 
degree 1. With this new function basis, we obtain a price of 4.876 for N = 60, 5 = 0.6 and 
M = 2 16 = 65536, within a computational time of 15 seconds. We note that this algorithm 
can be quite efficient in high dimension, even with the simplest choice of functions basis. 

Conclusion 

We have proposed a numerical scheme using empirical regressions on function bases. Re- 
garding the model, it is remarkably flexible and allows us to solve generalized BSDEs. 
Here the stochastic integral in the BSDE is driven by a Brownian motion, but we are 
preparing an extension of our scheme to deal with general martingales including jumps. 
Note that the terminal condition of the BSDE is here of the form but actually for 

certain path-dependences, our algorithm may still work after a slight modification (see 
[9]). Our approach is also suitable for reflected BSDEs. However, the robust choice of the 
penalization and regularization parameter as a function of the time step h is a delicate 
issue, which should deserve a special attention in future works. 
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